{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 变分量子线性求解器\n",
    "*Copyright (c) 2022 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*\n",
    "## 背景介绍\n",
    "\n",
    "线性方程组是数学中一个基本但非常有用的工具。 一个例子是，在经济学中，可以使用线性方程对经济进行建模。 此外，它还为非线性的大型系统提供了简单的估计。 因此求解线性方程组是一项重要的任务。\n",
    "\n",
    "变分量子线性求解器（Variational quantum linear solver, VQLS）是一种求解线性方程组的变分量子算法，采用了经典-量子混合的方案，可以在近期的含噪中等规模量子计算机上运行。具体来说，对于一个矩阵 $A$ 和一个向量 $\\boldsymbol{b}$，我们的目标是找到一个向量 $\\boldsymbol{x}$ 使得 $A \\boldsymbol{x} = \\boldsymbol{b}$. 使用 VQLS 算法可以得到一个与 $\\boldsymbol{x}$ 成比例的量子态，即一个归一化的向量 $|x\\rangle = \\frac{\\boldsymbol{x}}{\\lVert \\boldsymbol{x} \\rVert_2}$。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 模型原理\n",
    "\n",
    "量子场景的线性方程求解问题和通常的设定略有不同，因为量子计算需要将酉算子应用到量子态上。对于输入的矩阵 $A$，我们需要将其分解成酉算子的线性组合 $A = \\sum_n c_n A_n$，其中每个 $A_n$ 都是酉算子，可以在量子线路上运行。对于输入的向量 $\\boldsymbol{b}$，我们需要假设它是一个能够被某个酉算子 $U$ 制备的量子态 $|b\\rangle$，即 $U|0\\rangle = |b\\rangle$。我们可以用下面这张图来概括 VQLS 算法的整体架构:\n",
    "\n",
    "![VQLS](vqls.png)\n",
    "\n",
    "可以看到，VQLS 算法是一种混合优化算法，可以分为经典和量子两部分，需要在量子计算机上准备参数化量子电路 $V(\\alpha)$ 并计算损失函数 $C(\\alpha)$，然后在经典计算机上对参数 $\\alpha$ 进行优化从而最小化损失函数，直到损失低于某个阈值，最后输出目标量子态 $|x\\rangle$。其中参数化电路 $V(\\alpha)$ 可以生成一个量子态 $|\\psi(\\alpha)\\rangle$，电路 $F(A)$ 可以计算 $A|\\psi(\\alpha)\\rangle$ 与 $|b\\rangle$ 的近似程度，即损失函数 $C(\\alpha)$。当量子态 $A|\\psi(\\alpha)\\rangle$ 与 $|b\\rangle$ 足够接近时，这就意味着量子态 $|\\psi(\\alpha)\\rangle$ 与目标态 $|x\\rangle$ 足够接近，我们可以输出量子态 $|\\psi(\\alpha)\\rangle$ 作为目标态 $|x\\rangle$ 的近似。"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 量桨实现\n",
    "\n",
    "我们使用量桨中的 `Circuit` 类结合飞桨优化器来实现 VQLS 算法，其中量子部分中参数化量子电路 $V(\\alpha)$ 为 `Circuit` 中内置的 `complex_entangled_layer` 模板，损失函数计算电路 $F(A)$ 由 Hadamard Test 或 Hadamard-Overlap Test 组成，主要使用了量桨中的 `oracle` 量子门来实现控制 $A_n$ 门，在经典优化部分中我们使用 Adam 优化器来最小化损失函数。\n",
    "\n",
    "用户可以使用 toml 文件指定算法的输入，矩阵 $A$ 和向量 $\\boldsymbol{b}$，分别以 `.npy` 文件形式存储。用户可以使用以下代码，通过改变$n$的值随机生成一个 $n\\times n$ 的矩阵 $A$ 以及向量 $\\boldsymbol{b}$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "这是一个随机生成的A:\n",
      "[[4.1702199e+00+7.203245j   1.1437482e-03+3.0233257j\n",
      "  1.4675589e+00+0.9233859j  1.8626021e+00+3.4556072j\n",
      "  3.9676747e+00+5.3881674j ]\n",
      " [2.0445225e+00+8.781175j   2.7387592e-01+6.704675j\n",
      "  4.1730480e+00+5.5868983j  1.4038694e+00+1.9810148j\n",
      "  8.0074453e+00+9.682616j  ]\n",
      " [8.7638912e+00+8.946067j   8.5044211e-01+0.39054784j\n",
      "  1.6983042e+00+8.781425j   9.8346835e-01+4.2110763j\n",
      "  9.5788956e+00+5.3316526j ]\n",
      " [6.8650093e+00+8.346256j   1.8288277e-01+7.5014434j\n",
      "  9.8886108e+00+7.4816566j  2.8044400e+00+7.892793j\n",
      "  1.0322601e+00+4.4789352j ]\n",
      " [2.8777535e+00+1.3002857j  1.9366957e-01+6.7883554j\n",
      "  2.1162813e+00+2.6554666j  4.9157314e+00+0.5336254j\n",
      "  5.7411761e+00+1.4672858j ]]\n",
      "这是一个随机生成的b:\n",
      "[4.191945 +6.852195j  3.1342418+6.9232264j 6.9187713+3.1551564j\n",
      " 9.085955 +2.9361415j 5.8930554+6.9975834j]\n"
     ]
    }
   ],
   "source": [
    "n = 5\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "np.random.seed(1)\n",
    "A = np.zeros([n, n], dtype=\"complex64\")\n",
    "b = np.zeros(n, dtype=\"complex64\")\n",
    "for i in range(n):\n",
    "    for j in range(n):\n",
    "        x = np.random.rand() * 10\n",
    "        y = np.random.rand() * 10\n",
    "        A[i][j] = complex(x, y)\n",
    "    x = np.random.rand() * 10\n",
    "    y = np.random.rand() * 10\n",
    "    b[i] = complex(x, y)\n",
    "np.save(\"./A.npy\", A)\n",
    "np.save(\"./b.npy\", b)\n",
    "print(\"这是一个随机生成的A:\")\n",
    "print(A)\n",
    "print(\"这是一个随机生成的b:\")\n",
    "print(b)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "用户可以在 toml 文件中指定 VQLS 算法的参数 `depth`，`iterations`，`LR` 以及 `gamma`，分别对应参数化量子电路 $V(\\alpha)$ 的层数，优化器的迭代次数，优化器的学习率，和损失函数的阈值。在命令行输入 `python vqls.py --config config.toml` 即可完成线性方程组求解。这里我们给出一个在线演示的例子，首先定义配置文件的内容如下，用户可以自行更改 `test_toml` 中的参数："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_toml = r\"\"\"\n",
    "# 存储矩阵A的.npy文件的路径。\n",
    "A_dir = './A.npy'\n",
    "# 存储向量b的.npy文件的路径。\n",
    "b_dir = './b.npy'\n",
    "# 参数化量子电路的层数。\n",
    "depth = 4\n",
    "# 优化器迭代次数。\n",
    "iterations = 200\n",
    "# 优化器的学习率。\n",
    "LR = 0.1\n",
    "# 损失函数的阈值。默认为0。\n",
    "gamma = 0\n",
    "\"\"\""
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "运行 VQLS 算法如下："
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\yuzhan01\\Miniconda3\\envs\\pq_model\\lib\\site-packages\\paddle\\tensor\\creation.py:125: DeprecationWarning: `np.object` is a deprecated alias for the builtin `object`. To silence this warning, use `object` by itself. Doing this will not modify any behavior and is safe. \n",
      "Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations\n",
      "  if data.dtype == np.object:\n",
      " 88%|████████▊ | 176/200 [02:04<00:16,  1.42it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Threshold value gamma reached, ending optimization\n",
      "这是求解Ax=b的x: [ 1.3475237 -0.7860472j   0.22970617-0.88826376j -0.35111237-0.31225887j\n",
      "  0.07606918+1.2138402j  -0.729564  +0.48393282j]\n",
      "实际b的值: [4.191945 +6.852195j  3.1342418+6.9232264j 6.9187713+3.1551564j\n",
      " 9.085955 +2.9361415j 5.8930554+6.9975834j]\n",
      "算法得到的Ax的值: [4.185339 +6.8523855j 3.1297188+6.923625j  6.924285 +3.1467872j\n",
      " 9.092921 +2.932943j  5.8879805+6.999589j ]\n",
      "相对误差:  0.0008446976\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "import argparse\n",
    "import os\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\")\n",
    "os.environ[\"PROTOCOL_BUFFERS_PYTHON_IMPLEMENTATION\"] = \"python\"\n",
    "\n",
    "import toml\n",
    "import numpy as np\n",
    "import paddle\n",
    "from paddle_quantum.data_analysis.vqls import compute\n",
    "\n",
    "paddle.seed(0)\n",
    "\n",
    "if __name__ == \"__main__\":\n",
    "    config = toml.loads(test_toml)\n",
    "    A_dir = config.pop(\"A_dir\")\n",
    "    A = np.load(A_dir)\n",
    "    b_dir = config.pop(\"b_dir\")\n",
    "    b = np.load(b_dir)\n",
    "    result = compute(A, b, **config)\n",
    "\n",
    "    print(\"求解 Ax=b 的x:\", result)\n",
    "    print(\"实际 b 的值:\", b)\n",
    "    print(\"算法得到的 Ax 的值:\", np.matmul(A, result))\n",
    "    relative_error = np.linalg.norm(b - np.matmul(A, result)) / np.linalg.norm(b)\n",
    "    print(\"相对误差: \", relative_error)\n",
    "    np.save(\"./answer.npy\", result)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 引用信息\n",
    "\n",
    "```\n",
    "@misc{bravo-prieto2020variational,\n",
    "  title = {Variational {{Quantum Linear Solver}}},\n",
    "  author = {{Bravo-Prieto}, Carlos and LaRose, Ryan and Cerezo, M. and Subasi, Yigit and Cincio, Lukasz and Coles, Patrick J.},\n",
    "  year = {2020},\n",
    "  month = jun,\n",
    "  number = {arXiv:1909.05820},\n",
    "  eprint = {1909.05820},\n",
    "  eprinttype = {arxiv},\n",
    "  doi = {10.48550/arXiv.1909.05820}\n",
    "}\n",
    "```\n",
    "\n",
    "## 参考文献\n",
    "\n",
    "[1] “Variational Quantum Linear Solver: A Hybrid Algorithm for Linear Systems.” Carlos Bravo-Prieto, Ryan LaRose, Marco Cerezo, Yigit Subasi, Lukasz Cincio, Patrick J. Coles. arXiv:1909.05820, 2019."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pq-dev",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.15 (default, Nov 10 2022, 13:17:42) \n[Clang 14.0.6 ]"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "5fea01cac43c34394d065c23bb8c1e536fdb97a765a18633fd0c4eb359001810"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
